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! Abstract 

We apply the projection operator formalism to the problem of determining the 
asymptotic behavior of the lattice BGK equation in the hydrodynamic limit. As an 
alternative to the more usual Chapman-Enskog expansion, this approach offers many 
benefits. Most remarkably, it produces absolutely exact, though non-Markovian, hydro- 
dynamic difference equations as an intermediate step. These are accurate to all orders 
in Knudsen number and hence contain all of the physics of the Burnett equations and 
beyond. If appropriate, these equations may then be Taylor expanded to second order 
in Knudsen number to obtain the usual hydrodynamic equations that result from the 
J> ' Chapman-Enskog analysis. The method offers the potential to derive hydrodynamic 

■ difference equations for complex fluids with sharp gradients, such as immiscible and 
, amphiphilic flow, for which the assumptions underlying the Chapman-Enskog approach 

CN I are generally invalid. 

d ' 

^ ; 1 Introduction 

Statistical physics is often concerned with the problem of determining a closed set of equa- 
^ I tions of motion for a relatively small number of macroscopic degrees of freedom, given those 

■ for a much larger number of underlying degrees of freedom. This problem occurs in the 
classical theoretical context of deriving fluid dynamic equations from underlying kinetic 
equations [8], and also in the more modern computational context of so-called multiphysics 
simulations [131 [7]. These examples span a wide range of difficulty. 

The derivation of the Boltzmann equation assumes excellent separation between three 
different scales of length and time: 

1. The scales associated with a molecular collision are assumed to be the smallest and 
fastest of all the relevant scales. In particular, the range of intermolecular force is 
assumed to be much smaller than a mean-free path, and the duration of a collision is 
assumed to be much smaller than a mean-free time. 

2. The scales associated with the spatial and temporal intervals between collisions - the 
mean-free path and the mean-free time, respectively - are assumed to be much larger 
than the collision duration, but smaller than any hydrodynamic length and time scales. 
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3. Hydrodynamic length and time scales are assumed to be longest and slowest. 

Under these assumptions, an asymptotic approach may be adopted. This is the basis of 
the Chapman- Enskog expansion [8] with which it is possible to derive the Navier-Stokes 
equations of viscous hydrodynamics. 

Over the past few decades, it has been recognized that a breakdown of scale separation 
between scales 1 and 2 is much less catastrophic than a breakdown between scales 2 and 3. In 
a dense gas or liquid, for example, the mean-free path may be comparable to the interaction 
range, so good separation is lost between scales 1 and 2. This invalidates Boltzmann's 
stosszahlansatz, by which the Boltzmann equation is derived. In spite of this, as long as 
there remains good separation between scales 2 and 3, kinetic ring theory shows that the 
form of the Navier-Stokes equations is robust, and the interparticle correlations introduced 
by the loss of separation between scales 1 and 2 may be accounted for by an appropriate 
renormalization of the hydrodynamic transport coefficients [6]. The Navier-Stokes equations 
hold reasonably well for water at standard temperature and pressure, after all, even though 
scales 1 and 2 are comparable. 

A breakdown of separation between scales 2 and 3 is a much more serious issue. The ratio 
of mean-free path to hydrodynamic scale lengths is called the Knudsen number, Kn, and 
the Chapman- Enskog analysis is asymptotic in this quantity. Even for single-component, 
single-phase fluids, the Navier-Stokes equations may break down spectacularly when this 
quantity becomes too large. For complex fluid conflgurations, such as immiscible flow or 
two-phase coexistence, spatial gradients of order parameters may be very large indeed. The 
characteristic width of the interface between two immiscible fluids, for example, may be on 
the order of a mean-free path, so that the corresponding Knudsen number is of order unity. 
In this circumstance, asymptotic methods are not a viable option. 

For the last decade, physicists, chemists, applied mathematicians and engineers faced 
with the problem of modeling complex fluids have studied lattice models of hydrodynamics. 
These models consist of particle populations moving about on a lattice and colliding at 
lattice sites, whose emergent hydrodynamic behavior is that of a Navier-Stokes fluid [16]. It 
was found much easier to introduce effective forces between particle populations on a lattice 
than to introduce such forces in a continuum setting [15]. In this way, the dynamics of 
immiscible [9] , coexisting [1] , and amphiphilic [5] fluids have all been modeled successfully. 

It may be argued that the success of lattice models of fluids is purely phenomenological 
in nature. Attractive or repulsive interparticle potentials are introduced to make immiscible 
species separate, and surface tension is emergent. Potentials with both attractive and repul- 
sive regions, in the spirit of the Van der Waals potential, are introduced, and the liquid-gas 
coexistence is emergent. Simplifled BGK-style collision operators are used [2]. As long as 
the dimensionless parameters associated with the simulation match the fluid being modeled, 
the details of the microscopic interactions are deemed unimportant. 

Faith in this phenomenological approach is, to some extent, justifled by the observed 
robustness of hydrodynamic equations to details of the kinetic interactions. If both the 
interparticle potential range and the mean-free path are on the order of a lattice spacing, loss 
of separation between scales 1 and 2 is evident. As noted above, however, the hydrodynamic 
equations for a simple fluid are robust with respect to this loss; the hope is that a similar 
thing is true for more complex fluids. 
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Another reason for faith in the lattice-BGK approach owes to its second-order accuracy. 
Mathematically, the lattice-BGK equation may be written 



fj{r + cj,t + l)= fJ{r,t) + - 



^t^ ipir,t))-J){r,t) 



where fj(r,t) is the discrete- velocity distribution function corresponding to the j'th velocity 
Cj at spatial position r and time t, likewise /j"''^ is an equilibrium distribution function that 
depends only on the conserved densities p, and r is a collisional relaxation time. Second- 
order accuracy is not at all obvious from a cursory inspection of Eq. ([1]). To make it more 
evident, Dellar [10] has pointed out that we may define the new dependent variable 



(2) 



It is then straightforward to derive the lattice BGK equation for the transformed variable. 
Using P to denote the propagation operator, defined by 



Pfj{x,t) = f,{x + Cj,t+l] 



the result may be written 



T — 



I + P 



fr - fj 



(3) 



(4) 



where / is the identity operator. This form makes manifest the fact that the collision is 
applied between sites. It also allows a glimpse at the origin of the term r — 1/2 which will 
emerge as a factor in the transport coefficient. 

One might hope that the success of lattice BGK models would represent progress toward 
the end of deriving hydrodynamic equations for complex fluids from first principles. For 
example, it would be very satisfying if a Chapman-Enskog analysis of an interacting-particle 
lattice-BGK equation gave rise to a Ginzburg-Landau or Cahn-Hilliard equation for phase 
separation, coupled to the Navier-Stokes equations in the style of Halperin and Hohenberg's 
Model H [12]. To date, however, a successful derivation of this sort does not exist, most 
likely due to the aforementioned loss of separation between scales 2 and 3 for such fiuids. 
There is simply no small parameter analogous to the Knudsen number on which to base an 
asymptotic expansion. 

Part of the problem is that the very first step of the Chapman-Enskog analysis of the 
lattice-BGK equation is the Taylor expansion of the propagation operator, effectively in 
powers of the Knudsen number. This has the effect of yielding partial differential equations 
(PDEs) in space and time. In this work, we argue that the insistence on PDEs is the 
real problem. If we were willing to entertain the possibility of hydrodynamic equations 
that are difference equations in space and/or time, we would be able to model fluids with 
rapidly varying hydrodynamic flelds. The question arises: Is it possible to derive a closed- 
form equation for the hydrodynamic degrees of freedom of a lattice BGK model without 
Taylor expanding, or doing anything else tantamount to an asymptotic expansion in Knudsen 
number? In this paper, we answer this question affirmatively. 
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We circumvent the usual need for Taylor expansion by applying the projection operator 
formalism to the problem of deriving the exact hydrodynamics of the lattice-BGK equa- 
tion. As an alternative to the more usual Chapman-Enskog expansion, this approach offers 
many benefits. Most remarkably, it produces absolutely exact, though non-Markovian, hy- 
drodynamic difference equations as an intermediate step. These are accurate to all orders in 
Knudsen number and hence contain all of the physics of the Burnett equations and beyond 
for the lattice BGK fluid. If appropriate - but only if appropriate - these equations may then 
be Taylor expanded in space and/or time (to second order in Knudsen number) to obtain 
the hydrodynamic equations that would have resulted from the Chapman-Enskog analysis. 
The method thereby offers the potential to derive hydrodynamic difference equations for 
complex fluids with sharp gradients, such as immiscible and amphiphilic flow, for which the 
assumptions underlying the Chapman-Enskog approach are generally invalid. 

We begin by applying the methodology to a lattice-BGK model for Burgers' equation. 
This example is simple enough to display every step of the calculation, but complicated 
enough to raise most all of the above issues. In particular, spatially varying initial conditions 
may lead to a shock of characteristic width equal to the square root of the viscosity. For 
reasonable values of the viscosity, this shock width may be on the order of a lattice spacing. 



2 Burgers' equation 



2.1 Lattice BGK model 

The method is best illustrated by example, so we begin by applying it to a lattice BGK 
model for Burgers' equation in one spatial dimension [D = 1). A number of such models for 
Burgers' equation have been developed over the years [3], [11], 0]; this one is a variant of an 
entropic version considered by Boghosian et al. [1]. At each point x G Z, and at each time 
t G Z"*", we have a two-component distribution function f±{x,t), from which it is possible to 
recover the hydrodynamic density 

p{x,t) = U{x,t) + U{x,t). (5) 

The distribution function obeys the lattice BGK kinetic equation 

/±(x±l,t+l) = /±(x,t) ^ 



r L 



ft^ (p(x,t))-/±(x,t) 



where we have defined the local equilibrium distribution function 



(6) 
(7) 

That is, as the 



Here k is a parameter that is taken to be of first order in the scaling limit 
number of lattice points per physical distance is doubled, k will be halved. 

For the analysis of this model, we shall also find it useful to define the kinetic moment 

e(x,t):=/+(x,t)-/_(x,t), (8) 

so that the distribution function may be recovered from knowledge of its hydrodynamic and 
kinetic moments, 

p{x,t) ±^{x,t) 



f±ix,t) 



(9) 
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2.2 Chapman-Enskog analysis 

The simple lattice BGK model described above has pedagogical value as a simple introduction 
to the Chapman-Enskog asymptotic procedure. To make this presentation self-contained, 
and to facilitate comparison of the proposed new methodology with the Chapman-Enskog 
procedure, we outline the latter in this subsection. 

The usual analysis begins with a Taylor expansion of the kinetic equation ([6]) assuming 
parabolic ordering, wherein spatial derivatives are first order quantities, and time derivatives 
are second order quantities. The result is 

e'd,U ± edj^ + plh + O (e^) = ^ [f ± f P (l " " f±] ' (10) 

where all quantities are understood to be evaluated at (x, t), and where e is a formal expansion 
parameter, introduced for purposes of bookkeeping, that will be set to one at the end of the 
calculation. Note that we made the substitution k ento reflect the fact that k, is first-order 
in the scaling limit. 

Taking the sum of the ± components of this equation yields the conservation equation 

e%p + ed,^ + plp + O{e')=0. (11) 

To close this equation, it will be necessary to express in terms of p. This is done by solving 
Eq. (fTOj) perturbatively by taking 



f± = /f + + e'fi'^ + O {e') . (12) 
At order zero in e, we immediately obtain 

/f = f . (13) 

This lowest approximation to f± implies = 0; it is therefore insufficient to calculate the 
kinetic moment, which first appears at order one. 
Proceeding to order one, we obtain 

±a./f = i[±^p(i-0-/i"], (14) 
/^'^±ip(i-f).^a... (15) 

where we have used the order zero solution in the last step. This order-one result for the 
distribution function, 

p en / p\ €T 

yields the kinetic moment 



/.^i±fp(l-ij^^^.P + 0(.^) (16) 



^ = e^p{l-^)-eTd,p+0{e'). (17) 
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Inserting this result in Eq. ffTTl) yields 



e^dtp + ed^ enp ( 1 



P 
2 



) 



p + O (e=^) = 



(18) 



or 



dtp + K (1 - p) d^p 




p + 0{e). 



(19) 



Upon the substitution u := k (1 — p), this becomes Burgers' equation in canonical form, 



where we have defined the transport coefficient v := t — 1/2. 

This method of simulating Burgers' equation is simple to implement and remarkably 
robust. Fig. [1] shows the results of simulating of the model on a periodic lattice of size 
N = 1024, with K, = 0.25 and r = 0.60. The initial conditions used were 



It is seen that the solution captures the formation and decay of the shock, and that the 
width of the shock at late times is comparable to the lattice spacing. 

Before going on to the exact analysis, it is worth making some general observations about 
the Chapman-Enskog analysis for this model: 

• We solved the kinetic equation only to first order, but used that solution in the hydro- 
dynamic equation at second order. This interlacing of orders is characteristic of the 
Chapman-Enskog analysis. 

• The second-order accuracy of the lattice BGK equation is not at all evident in the 
final result obtained. If we were to carry on to the next order - which would involve 
solving the kinetic equation to second order - it is not at all clear that we would not 
find corrections to the hydrodynamic equation obtained. 

• The transport coefficient is equal to r — 1/2. The first term, r, arose from the gradient 
correction to the local equilibrium distribution function. The second term, —1/2, came 
from the Taylor expansion of the left-hand side of the kinetic equation at order two, 
and is an artifact of the discrete nature of the model. Note that the combination 
r — 1/2 is manifest in Eq. (j4]). 

2.3 Exact analysis 

The exact analysis that is the point of this paper comes from projecting the kinetic equation 
onto hydrodynamic and kinetic subspaces, solving for the kinetic field ^ (x, t) as though the 
hydrodynamic field p(x, t) were a known function of position and time, and then using this 
solution to obtain a closed dynamical equation for p{x,t). This general technique has been 
known for a long time under a variety of different names. In physics, it is sometimes called 
the Mori-Zwanzig projection operator formalism [19]. In linear algebra, it is related to the 



dtu + u dxU = V d\u + C (e) , 



(20) 




(21) 
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Figure 1: Simulation of lattice BGK model for Burgers' equation. The formation and subse- 
quent decay of the shock are evident. Note that the shock width is on the order of a lattice 
spacing, calling into question the very first step of the Chapman-Enskog analysis, namely 
the assumption that spatial gradients are small. 
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Schur complement [18]. To the best of our knowledge, however, it has not heretofore been 
applied to the problem of obtaining exact solutions of the lattice BGK equation. 

We begin by using Eqs. ( P[6|9|) to write coupled evolution equations for p and ^. After a 
bit of straightforward algebra, we obtain 



+^[p{x + l,t)+p{x-l,t)] 



+ 27 



-p{x + l,t) 1 



2r 

p(x + l,t) 



[^{x + l,t)-^{x-l,t)] 



l [p{x + 1, t) - p{x - 1, t)] + ^-^ [^{x + 1, t) + e(x - 1, t)] 
z Zt 



K 

+ 27 



+p{x + l,t) 1 



p(x + l,t) 



+ p(x-l,t) 1 



p{x - 



(22) 



.(23) 



These two dynamical equations for the hydrodynamic and kinetic fields respectively, taken 
together, are equivalent to the original kinetic equation. The equation for the kinetic field 
may be written in the suggestive form 



i{x,t + l] 



2r 



[e(x + l,t)+e(x-l,t)] + Q(x,t), 



(24) 



where we have grouped together everything involving the hydrodynamic field into a single 
effective source term, defined by 



[p(x + l,t)-p(x-l,t)] 



p(x + l,t) 1 



p{x + l,t) 



+ p(x-l,t) 1 



p{x - l,t) 



(25) 



Eq. 



is a linear difference equation. For an infinite lattice, it has the exact solution 

t-\ n 



e(M)^EE(::)(i)"« 

n=0 m=0 ^ ^ 



X + 2m — n,t — n) 



+ (|)*e(a: + 2m-t,0), (26) 

m=0 ^ 



where we have defined a := 1 — 1/r, and we note that |(t| < 1 for r > 1/2. The first term 
above describes the excitation of the kinetic field due to gradients of the hydrodynamic field. 
The second term is a transient, due to the initial conditions. 



2.4 Self-consistent hydrodynamic difference equation 

Supposing, for the time being, that the kinetic field is initialized to zero 0, or that we have 
waited long enough for transient behavior to be unimportant, we can insert the first term 
on the right-hand side of Eq. (l26l) into the dynamical equation for p and rearrange to obtain 

^It should be evident that this is an inessential assumption that we make here only for the sake of 
simplicity of presentation. 
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the nonlinear difference equation 



1 



p{x, t + 1)- p{x, t) = +- [p{x + 1, t) - 2p(x, t) + p{x - 1, t)] 



2r 



p(x + l,t) 



2 



p{x - l,t) 



^ _ p{x - l,t) 
2 



2, 



n=0 'm=0 

[p{x + 2m — n + 2, t — n) — 2p(x + 2m — n,t — n) + p{x + 2m — n — 2,t — n)] 

t-l n 



p{x + 2m — + 2, t — ra) 



n=0 r?i=0 

-p(x + 2m — n — 2,t — n) 



p{x + 2m — n + 2, t — 



p{x + 2m — n — 2, t — n) 



(27) 



The remarkable thing about this result is that no essential approximations have been made 
in its derivation. It is an exact difference equation that must be obeyed by a hydrodynamic 
field p(x, t) satisfying the lattice BGK equation. 



2.5 Recovery of a hydrodynamic differential equation 

Examination of Fig. [1] indicates that the hydrodynamic field changes reasonably slowly in 
time so that a Taylor expansion in time is justified, but that it can change very rapidly 
in space after the onset of the shock, rendering a spatial Taylor expansion of questionable 
validity. For this reason, we begin by Taylor expanding in time only, incurring an error of 
at most second order, 



dtp{x, t) = +^ [p(x + 1, t) - 2p(x, t) + p(x - 1, t)] 



p{x + l,t) 



p{x - l,t) 



p{x - 



t~l n 



m) \2. 



n=0 m=0 

[p(x + 2m — 77, + 2, t) — 2p(a; + 2m — n, t) + p{x + 2m — n — 2, t)] 

t-l n 



K 



EE 



n=0 m=0 

-p(x + 2m — n — 2, t) 



p{x + 2m — n + 2, t) 



1 - 



p(x + 2m — n — 2,t) 



p{x + 2m — n + 2, t) 
2 

- o (^) ■ 



(28) 



This is a second-order accurate set of coupled discrete-space, continuous-time equations 
describing the hydrodynamics of the model. 
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To proceed to a spatiotemporal partial differential equation, we may Taylor expand 
Eq. fl28l) about the spatial point x, retaining spatial derivatives to second order. Because 
of the symmetry of the differences, it is straightforward to see that the error incurred is of 
second order. The sum over m is then a binomial series and that over n is a finite geometric 
series, yielding 



dip + O (e^) . (29) 



Ignoring the terms that decay exponentially in time (note |cr| < 1), we recover Eq. f|T9|) . 

Before going on to consider the exact analysis for more general lattice BGK equations, 
we offer some preliminary observations: 

• The usual Chapman-Enskog procedure begins with a Taylor series expansion of the ki- 
netic equation in difference form, and then assembles continuum hydrodynamic equa- 
tions from that series. This method, by contrast, first derives exact hydrodynamic 
equations in difference form, and only then (optionally) Taylor expands the results to 
obtain differential equations. 

• Closely related to the above point is the idea that the hydrodynamic difference equation 
first obtained is not expanded in Knudsen number, and therefore expected to hold 
in limits that would not ordinarily be considered "hydrodynamic." Such limits may 
include situations with steep spatial gradients or long mean-free paths. 

• The method yields an exact hydrodynamic equation in difference form that can be Tay- 
lor expanded to yield the hydrodynamic equation. The second-order accurate nature 
of this expansion is manifest. 

• The method works only for lattice BGK equations, and relies on the fact that the 
equilibrium distribution function is a function of the conserved quantities only. (If this 
were not the case, Eq. (1211) would not be linear in the kinetic field, so we would not be 
able to solve it exactly.) 

• Note that as t ^ oo, the effect of the last two terms on the right-hand side of Eq. fl28l) 
is to simply alter the coefficients in front of - or "renormalize" - the other terms in 
the equation. The penultimate term evaluates to Cpxx so it renormalizes the diffusion 
term, and the last term is —Dnpil — p)px so it renormalizes the advection term. The 
sum determines that C = ar and D = a whence 

9tP = ^Pxx - np (l - + CT'rpxx - crdx up (l - , (30) 
from which Eq. (fT9l) follows. 



3 General procedure 
3.1 Projection operators 

We now suppose that we have a 6-component discrete- velocity distribution function fj{r,t), 
where j £ {1, . . . ,b}, where r G £ is a point on lattice C, and where t G Z"*" is the time. The 
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bh hydro dynamic variables are obtained by projection with the x b matrix H, 



p„ = Y.^c,'fv (31) 
j 

and the bk = b — bh kinetic variables are obtained by projection with the bk x b matrix K, 

b 

= EV/.- (32) 

j 

The rows of H and K are linearly independent, so that knowledge of all the hydrodynamic 
and kinetic variables is sufficient to reconstruct the distribution function, 

/.• = E^.% + E^A- (33) 
It is easily verified that the b x bh matrix A and the b x bk matrix B obey the relations 

b b 

E^"V = ^/ E^«'^/ = (34) 
j j 

EVV = EV'^/ = V' (35) 

j j 

and 

V = E^/V + E^/^.'- (36) 

In the above presentation, we have adopted the convention of using Greek letters from the 
beginning of the alphabet (a, /3, . . .) to label the hydrodynamic variables, from the middle of 
the alphabet (yU., i', . . .) to label the kinetic variables, and Latin letters {j, k, . . .) to label the 
distribution function components. We shall adhere to this notational convention as closely 
as possible in the forthcoming development. 

Example 1 A triangular grid in two spatial dimensions (D = 2) has the 6 = 6 lattice vectors 

Cj := cos I — 1 + By sm I — 1 . (37) 

In the widely adopted nomenclature for lattice Boltzmann models, this is called the D2Q6 
model [§. // mass and momentum are conserved, the hydrodynamic and kinetic projection 



^We note in passing that the D2Q6 model is no longer extensively used for two-dimensional lattice BGK 
simulations of fluids. It has been abandoned in favor of the so-called D2Q9 lattice which may be implemented 
on a Cartesian grid. We use the D2Q6 model in this example only for the sake of simplicity, since it allows 
us to exhibit matrices of dimension six rather than nine. It should be clear that there is nothing preventing 
this method from being applied to the D2Q9 lattice, or even to the much larger lattices used in modern 
lattice-BGK simulations. 



11 



operators may be taken to be 



H 



and 



K 



1 
1 


1 
1 



1 


1 


1 


1 


1 " 


1 


1 


-1 


1 


1 


2 


2 


2 


2 


V3 


V3 





V3 


V3 


2 


2 


2 


2 J 



(38) 



-1 


1 


-1 


1 


-1 " 


1 


1 


1 


1 


1 


2 


2 


2 


2 













2 


2 


2 


2 - 



(39) 



respectively. Note that the first row of H corresponds to the mass density p, while the second 
and third rows correspond to the momentum density tt. Collectively, we refer to the conserved 
densities as p = (p, tt) . The corresponding reconstruction matrices are then 



A 





6 
6 



6 

V3 



B 





V3 
6 

V3 
6 



6 

V3 
6 



(40) 



The identities of Eqs. ^^3^ through ^3^) are then readily verified. 



3.2 General form of the lattice BGK equation 

The general form of the lattice BGK equation, Eq. ([1]), may be rewritten 



/, (r, t + l)= af,{r - c„ t) + (1 - a) /J-^ (p (r - c„ t)) 



(41) 



where cj := 1 — 1/r is defined as before. Note that the equilibrium distribution function /^'"'^ 
is allowed to depend only on the hydrodynamic moments p. In what follows, it shall prove 
useful to write the equilibrium distribution in the form 



(42) 



Comparing this with Eq. 0331) . we see that Q (p) is the kinetic portion of the equilibrium 
distribution function. 

Example 2 For a fluid with p = (p, tt), the form generally used for this dependence is the 
Mach- expanded distribution 



f^''\p) = W, 



P + — TT ■ Cj + 



1 



77 TT 



: {cjt 



c%) 



(43) 
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where the Wj are weights associated with each direction, I2 is the rank-two unit tensor and 
Cs is the sound speed defined by the isotropy requirement 



(44) 



In particular, for the D2Q6 lattice of Example Ui it may be verified that Wj = 1 and Cg = 
1 / \/2. The first two terms of Eq. are then the hydrodynamic portion of the equilibrium 
distribution function, and 



(45) 



is the kinetic portion. 



By means of the projection operators defined in tlie previous subsection, it is straightfor- 
ward to decompose this into coupled evolution equations for the hydrodynamic and kinetic 
moments, 



Po.{r,t+l) = (tJ^H, 



and 



£ A/p^{r - c„ t) + Yl B.'Ur - c„t) 



e/.(r,t+l) = a 



b 

(l-^)E^V/)"^ (p(r-c„t)). 



(46) 



(47) 



respectively. Eqs. ( 146|) and (1471) are the generalizations of Eqs. ( 122|) and (l23l) . respectively. 
The equation for the kinetic modes may be written more succinctly as 



bk b 



(4J 



where we have defined 

bh b 



^) ^= E E (^^/'^/) p^^"- - ^) + (1 - ^) E ^1"^ ('^ - ^)) • (49) 

Eqs. ( l48l) and (l49l) are the generalizations of Eqs. (j24l) and (j25l) . respectively. 
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3.3 Exact analysis in the general case 

To proceed as in the example, it is now necessary to find an exact solution to the linear 
equation Eq. ( 148|) for the kinetic modes, assuming that the hydrodynamic modes are known. 

Consider a labeled path of n steps along lattice vectors, whose vertices are labeled by 
kinetic modes. More specifically, consider a path along the lattice beginning at position r' 
and mode yU.„ = at time t' = t — n, and ending at position r & C and mode /iq = /i at time 
t. One such path in the D2Q6 model is illustrated in Fig. [2l The set of all such n-step paths 
will be denoted by p(")(r', v; r, yu). 

A path p G V^"'\r' , u] r , fx) is thus characterized by its sequence of indices j{p) : = 
{jn, ■ ■ ■ ,ji} of the lattice vectors traversed, and also the sequence of kinetic modes /x : = 
{fin, . . . , f^o} at the visited vertices. Note that it must be true that 

n 

Y^c,^=r-r', (50) 



i=i 



and 



/io = At (51) 
/i„ = u. (52) 

We use Sp ' '^'''"''^^ to denote the sum over all paths, p G V^^^{r', u; r, //). 

Example 3 A path of length n = 5 in the D2Q6 model of Example [3 is shown in Fig. [3 
The sequence of lattice vector indices pictured is j{p) = {3, 2, 5, 1, 3}. 

To a path p we assign the weight 

n 

^^"^;(j,A^)=n(^..-/'^^./0' ^^^^ 

where there is no understood summation on repeated indices. The exact solution to Eq. (l48l) 
is then 

t-l bk V'-"'>(r',u;r,n) 

n=0 r'eC V p 

+ E E E ^' A^(P)) Ur', 0). (54) 

r'e£ u p 

Note that if r' can not be connected to r by a sequence of n lattice vectors, perhaps because 
it is too far away, then V^^\r', /i^; r, /xq) is understood to be the null set. 

Example 4 Eq. ( [5^ is the generalization of Eq. [2D\] . To see this, note that our D = 1 
example for Burgers' equation had bh = bk = 1 and 6 = 2. The projection operators were 

H = [+1 +1] (55) 
i^' = [ +1 -1 ] , (56) 
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Position r 




and 



from which it follows that w^'^\j . 



A 



B 



-1/2 
-1/2 

-1/2 
-1/2 



(57) 



(58) 



The number of paths connecting r' = x+2m- 



-n 



and r = X is then the binomial coefficient (^), resulting in Eq. ( f^) . 



3.4 Self-consistent hydrodynamic difference equations 

As with our one- dimensional example, we suppose that the kinetic field is initialized to 
zero, and we insert Eq. (!54l) into Eq. (146|1 and rearrange to obtain the exact hydrodynamic 
difference equation 

Pa{r,t + 1) 

b bh b 

= ^ E E ^a^'^/pMr- - c„t) + {l-a)Y. HJ ) (p (r - c„ t)) 

j /3 j 
b bk t-l r'-"'>(r',u;r-Cj,fi) 

+ E g.(r',t-n).(59) 

Eq. (1591) is an exact nonlinear hydrodynamic difference equation for the conserved den- 
sities, albeit in terms of a diagrammatic summation. We are now free to Taylor expand in 
either time or space, as appropriate to the phenomenon under consideration. Taylor expan- 
sion of Eq. 0591) to second order in the time variable followed by a bit of rearranging yields 
the continuum-time, discrete-space hydrodynamic equation 

dtpair,t) 

b bh 

= EE V ip^^"- - ^) - /^/^(^' ^)] 

b 

j 

b bk t-l Vf"\r',u;r-Cj,fi) 

+EEEE E (60) 

j fi,u n=0 r'£C p 

3.5 Recovery of a hydrodynamic differential equation 

Finally, to proceed to a hydrodynamic differential equation, we may Taylor expand in space, 
retaining terms to second order. While this step is not technically difficult, it is tedious 
enough to warrant splitting the calculation into three parts, corresponding to each of the 
three terms on the right-hand side of Eq. fl60l) . We begin by defining some useful tensors in 
terms of the projection and reconstruction operators and the lattice vectors. 



(P(r-c„t))-er iPir,t)) 
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3.5.1 Some useful tensors 



In preparation for the forthcoming analysis, it is useful to define the tensors 



N 



SfiN) := E^-'V®'^. (61) 



TJ{N) := E^^'ll®^^- (62) 



N 



US'(N) :^ J:hjJ^_^.„ (63) 

where ®^ denotes an A'"-fold outer product. If all of the conserved quantities are scalars, 
these tensors have spatial indices ranging from 1 to D in which they are completely 
symmetric; that is, all have spatial rank A^. In addition, the S{N) and T{N) each have two 
hydrodynamic indices ranging, a and /3, ranging from 1 to bh, and the U{N) have three such 
hydrodynamic indices. 

When it becomes necessary to refer to these tensors by their spatial indices, we adopt the 
convention of replacing the number in parentheses by the actual list of spatial indices. 
Thus, for example, we have 

b 

Sj{{a, b, c}) = E HJA/c,aC,bC,c, (64) 
j 

where Cja is the ath spatial component of lattice vector Cj. Note that the value of A^, namely 
A^ = 3 in this case, can be inferred by the fact that there are three indices in the list. When 
listing components for A^ = 0, we denote the empty list by {}. 

If one or more of the conserved quantities are spatial vectors, such as momentum, it is 
convenient to abuse notation by allowing the hydrodynamic indices to be either scalars or 
vectors. Each one that is a vector will increase the spatial rank of the tensors by one. For 
example, if ( is the index of a conserved scalar and rj is the index of a conserved vector, 
S^'^{N) is a completely symmetric tensor of spatial rank A^, while S^'^{N) is a tensor of 
spatial rank A^ + 1 that is symmetric under interchange of A^ of its indices, and S^'^{N) is a 
tensor of spatial rank A^ + 2 that is also symmetric under interchange of A^ of its indices. 

To specify components of these, it may be necessary to use nested subscripts. Thus, we 
have that 

U,J'({b^^}) = T.H,j£^^c,,c,. (65) 

is the (a, b, c) component of a tensor of spatial rank three that is symmetric under interchange 
of its rightmost two indices. Note that the notation automatically ensures that the A^ 
symmetric indices will be the rightmost indices. 

Finally, we also define alternative versions of the S{N), T{N) and U{N) tensors, using 
the same names for economy of notation, but with upper hydrodynamic and lower kinetic 
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indices, 

b N 

SfiN) := E^/^/®^^- (66) 

j 



U^'(N) :^ (68) 

All of the above-mentioned considerations about counting independent components likewise 
apply to these alternative versions. 

Example 5 For the D2Q6 model in Example [21 supposing that mass and momentum are 
both conserved, the matrices H , K, A and B are given in Eqs. / T^gj) through (^^. // we 
denote^ the hydrodynamic indices by a,P,'j E {p,7z}, it is straightforward to compute the 
above tensors. Tabulations of the results for S{N), T{N) and U{N) are shown in TablesUl 
m and\^ respectively. We have not bothered listing entries that are anisotropic and/or not 
needed for the derivation of the hydrodynamic equations. 

To assemble the hydrodynamic equations from Eq. (!60|) . we label the terms on its right- 
hand side as 0, @ and (3), and consider each separately. 

3.5.2 Evaluation of first term in Eq. (I60l) 

The first term on the right of Eq. fl60l) can be written 



®« ■= J2Y.^c.'A/[pp{r-c,,t)-pp{r,t)] 

j P 

b bh 



j 13 

bh 1 bh 



Cj ■ Vpfs{r,t) + ^CjCj : Wpf}{r,t) 



= -^5/(1) ■Vp^(r,t) + 1^^/(2) :VVp^(r,t) + --- (69) 

f3 (3 

In Appendix [XJ we evaluate this expression for the D2Q6 model. When doing so, we adopt 
incompressible scaling [T3], wherein spatial gradients and Mach number are taken to be order 
e, and time derivatives and density fluctuations are taken to be order e^. To leading order, 
we find that the p and tt components of the above term are 

®p = -V-TT (70) 

and 



■^Unfortunately, denoting mass by index p and momentum by vector index tt violates our convention 
that indices representing conserved quantities should come from the beginning of the Greek alphabet. The 
standard use of these letters to represent mass and momentum density, however, was deemed to override 
this concern. 
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(71) 



3.5.3 Evaluation of second term in Eq. (I60p 

Similar considerations can be applied to the second term on the right of Eq. fl60l) which may 
be written 



(72) 



or 



1 - a 



bh 



J2^Jil)-^Pp+lf2^Ji^) ■■ VVp.+lf^U/'i^) ■■ (Vp,)(Vp,)+--- (73) 



/3 13 /3,7 



In Appendix |Bl we evaluate this for the D2Q6 model in the limit of incompressible scal- 
ing [H] . To leading order, we find that the p and tt components of the above term are 



and 







^ f ^ 1^1 ,2 
TT ■ VtT V 77 

2p\ 2 ' ' 



(74) 
(75) 



3.5.4 Evaluation of third term in Eq. (I60p 

The third term on the right-hand side of Eq. (l60l) involves Q^{r, t) which is given by Eq. ( H9l) . 
so it is necessary to compute this first. We note that it may be expanded to second-order 
accuracy to obtain 



+ (1 - ^) 



P/3 



4Ef^/"(2)aVp,) (Vp,) + --- 



/3,7 



(76) 



In Appendix [CI we evaluate this for the D2Q6 model in the limit of incompressible scal- 
ing [13]. To second order, we find that its p and tt components are 



Qn 



TT - -Vp + -V^ 



77 ■ V77 V 77 



(77) 
(78) 
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As noted for the example of Burgers' equation, the effect of the diagrammatic sum as 
t — > cxD and in the continuum hmit is to apply a linear operator to these terms, thereby 
renormalizing other terms in the equation. The most general form for the third term in 
Eq. ([60]) is then 



Qa Ca 



1 



Vp+ -V^TT 



1 - a 



TT ■ VtT V IttI^ 



+ D^Vp + ^« W, (79) 



2 ' 8 \ 2p J \ 2 

where C^, and are determined by the diagrammatic series. The results are 



a = = E,, = 



and 



C71 



l + a 
1 - a 
1 fl + a 



2 V 1 - a 
1 / (T 



4\l-a 

so the third term of Eq. (!60|) has no p component, 

@p = 

and TT component equal to 



1 


+ 


a 


1 




a 




1 


(\ 


+ 


2 





--Vp+ -V^TT 



Vp+ - 



1 - a 



2p 



TT ■ VtT V IttI 



a 



4 V 1 - (T 



3.5.5 Assembly of hydrodynamic equations 

Armed with these terms, we may now assemble the full hydrodynamic equations, 

Pt = ©p + (2), + @, 
= ©. + (2). + ®.. 

Since pt may be ignored in the incompressible limit, the first of these reduces to 

V • 77 = 0, 

while the second gives 

l + a 



dtTZ 



1 



+ - 



1 - a 
1 /l + a 



2 8 



2p 



2 V 1 - (7 



Vp+ - 



a 



4 V 1 - a 
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(80) 

(81) 
(82) 

(83) 
(84) 



^5) 



(86) 
(87) 



^9) 



This simplifies to yield 

dt-K + -77 ■ VtT = - VP + Z/V^TT, (90) 

P 

where the pressure P is given by the equation of state 

P=^-^, (91) 

and the kinematic viscosity is 

1 f 1 



Eqs. (1551) and flMIl) are seen to be the Navier-Stokes equations of viscous, incompressible 
hydrodynamics. The expressions for the equation of state Band the viscosity are well known 
for the D2Q6 fluid 



4 Discussion 

Diagrammatic methods often yield new physical insights, and this one is no exception. Since 
the Qa are the driving terms in the linear equations for the kinetic modes we may 
think of the Qa as the precise combination of hydrodynamic gradients that excite kinetic 
modes. Thus excited, the kinetic modes may propagate and couple with other kinetic modes. 
The diagrams trace the evolution of these kinetic excitations in space and time until they 
ultimately project, via the B^^ matrices, a contribution back to the hydrodynamic modes. 
Alternatively stated, the sum over diagrams gives the Green's function of Eq. fHHj) which 
governs the evolution of the kinetic modes. 

It is remarkable that the effect of these kinetic excitations is often nothing more than the 
renormalization of terms already present in the hydrodynamic equation. The diagrammatic 
sum is necessary to obtain the correct answer for the transport coefficients (advection, diffu- 
sion and viscosity in the above examples), but it is not necessary to obtain the general form 
of the hydrodynamic equation. In addition, principles of symmetry and covariance may be 
employed to reduce the diagrammatic sum to the determination of a few scalar quantities - 
the A, B and C quantities in the above examples. 



5 Conclusions 

We have described a new method to derive hydrodynamic equations for lattice BGK fluids. It 
is qualitatively different from the usual Chapman-Enskog analysis, and superior insofar as it 

^The equation of state includes a term that is clearly nonphysical because it depends on the hydrodynamic 
velocity. Great concern is often expressed about this term, but it is entirely misplaced and indeed reflects 
a serious misunderstanding of the nature of the incompressible limit, in which the form of the equation 
of state is entirely irrelevant and where the pressure is given by the solution to a Poisson equation. As 
long as one takes care to keep the Mach number small enough so that the density fluctuation scales as its 
square, the method is perfectly valid for simulating incompressible flow, the form of the equation of state 
notwithstanding. 
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results in absolutely exact hydrodynamic equations. Additionally, while more demanding in 
terms of calculation, the method is arguably conceptually simpler than the Chapman- Enskog 
approach, easier to automate using symbolic mathematics software, and more transparently 
a description of the generation and propagation of kinetic modes in a moving fluid. 

We have illustrated this new methodology, first by presenting a simple lattice-BGK model 
for Burgers' equation, and second by presenting a more complex lattice-BGK model for a vis- 
cous, incompressible fluid. The final step in this method is the extraction of a diagrammatic 
sum, but we have shown that all that we really need is the asymptotic form of this sum for 
n large, and also that the effect of this sum is often nothing more than the renormalization 
of hydrodynamic transport coefficients. 

Because the method relegates the Taylor expansion of the propagation operator to an 
optional step at the end of the analysis, it is possible to derive hydrodynamic equations that 
are either discrete or continuous in space and/or time. For the examples presented in this 
work, we showed the result of expanding in time but not space. For certain phenomena in 
transport theory that tend to be uniform in space but rapidly changing in time, such as 
spontaneous micelleization, it may be more appropriate to expand in space but not time. 

Future work will take up the application of this method to complex fluid phenomena, such 
as multiphase fluids described by the Shan-Chen model [15]. It is hoped that this method 
will yield accurate hydrodynamic equations in the spirit of Halperin and Hohenberg's Model 
H [12] for such complex fluids. 
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A Appendix: Evaluation of Eq. ( 1691 ) for D2Q6 fluid 



For a mass and momentum-conserving ffuid, we wish to compute the p and tt components 
of Eq. fl69l) . The p component may be written 

b bh 

®p ■= Y.Y.Hp'Af[pp{r-c,,t)-pp{r,t)] 

+ ^S/{{a, 6}) V„V,p + ^S/^{{b, c})V,V,7r, + ■ ■ • (93) 
Using the resuhs for the D2Q6 model compiled in Tables [H [2] and [31 this becomes 



1 9 

= -V-7r + -V> + --- 

4 



(94) 



When incompressible scaling is taken into account, the second term is order times the first 
one, so we may write the result. 



-V ■ 77 + C (e^) . 



(95) 



The TT component may be written 

b bh 



= -^^/({6})V,p-S^/H{c})Ve7r, 



(96) 



Using the results for the D2Q6 model compiled in Tables [H [2] and [3], this becomes 



is^/iib, C}) V^VeP + Is^J'iic 4) VcV^TT, + 



/ o 



(97) 



When incompressible scaling is taken into account, the second term is order times the first 
one, so we may write the result. 



®. = Vp + i [W + 2V (V . tt)] + O (6^) . 



(9^ 
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B Appendix: Evaluation of Eq. (|73l) for D2Q6 fluid 



For a mass and momentum-conserving fluid, we wish to compute the p and tt components 
of Eq. (1731) . The p component may be written 



l-a 



Y.HAit\p{r-c,,t))-i^p\p{r,t)) 



-T/({a})V«p-T/"(W)V,7r„ 

+ ^T/({a, h})VaV,P + ^T/" ({6, c}) V.VeVr, 

+\u/P{{a,h}) (V.p) (V,p) + ^f//-({6,c}) (V,p) (VeTT.) 



p^ ^ ' ■ 4p 

^ (TT^^afc - T^aT^b) C^aP) (Vfep) - ^ (^afcVTc + ^acTTfe - SbcT^a) (Vfep) (VcVTa) 
1 



4p" 



(^afeVTc + ^acTTfe - ^fccTTa) (VfeVTa) (VcP) 



+ 1" (^ac^bd + Sad^bc - Sabred) cT^a) (VrfTTfe) + 

4p 



(99) 



Using the resuhs for the D2Q6 model compiled in Tables [H [2] and [3], this becomes 

|7r| VV - 2-7777 : VVp) + — [tt ■ V ( V ■ tt) + tt ■ V ( V ■ tt) - tt ■ V^tt] 



1 - (T 



5p- 



4p 



4p3 

1 

v 

1 



[|7rnVpr-(7r.Vpf] 



[(tt ■ Vtt) ■ (Vp) + (tt ■ Vp) (V ■ tt) - (Vp) ■ (Vtt) ■ tt] 
(tt ■ Vp) (V ■ tt) + (tt ■ Vtt) ■ (Vp) - ^ (V IttI^) ■ (Vp) 



+— {(V ■ tt)' + Tr [(Vtt) ■ (Vtt)] - (Vtt) : (Vtt)} + ■ ■ • 

All of these terms are negligible when incompressible scaling is applied, so we have 



(100) 



+ O (e^) . 



(101) 
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1 - a 



The TT component may be written 

/ON ^ 

- [^f''\pi.r - c„t)) - S}r\p{r,t)) 

j 

= -T^/({6})V,p-r^/H{c})Ve7r, 

+^[/^/^({6,c}) (V,p) (V,p) + if/^/-H{c,rf}) (V,p) (V^TT,) 
+ ^[/,/^^({c,rf}) (VeTT,) (V,p) + ^f/./^^1{rf,e}) (V,vr,) (V.vr,) + 
Using the resuhs for the D2Q6 model compiled in Tables [H [2] and [31 this becomes 

1 - (T 



2p 



■(102) 



(103) 



Applying incompressible scaling, we are left with the dominant terms 



1 - a 



TT • VtT + TtV ■ TT V IttI 



(104) 



C Appendix: Evaluation of Qa for D2Q6 fluid 

For a mass and momentum-conserving fluid, we wish to compute the p and tt components 
of Qa, given in Eq. ( 1761) . The p component may be written 

Qp = + V(W)VaP- 

+ ls/{{a, b})VaVkP + ^ C}) ViVeTT^ 

- (1 - a) T/({a}) V.p - (1 - a) T/"({6}) V.vr, 

T/({a, 6}) V.Vbp + i^T/''({6, c})V,VeVr, 



1 

H — 




cr 




2 




1 

+ - 




cr 




2 




1 




a 


' 2 



2 



(105) 
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Using the results for the D2Q6 model compiled in Tables [H [2] and [3], this becomes 

Qp = +P- ^ab^bT^a + \5ab^aSbP + {^"^^ab " '^T^a'^b) ^ bP 

4 8p^ 

1 (T 1 
H ^ (^afoVTc + ^acTTft - ^focVTa) VfoVcVT^ {rc'^^ab - T^aT^b) ("^ ap) (Vfep) 

4p2~ (^'^^^'^ + '^"^^^ ~ ^bcTTa) (Vfep) (VcVTa) 

-TIT- (^abVTc + (JacVTb - SbcHa) (VfeTTa) (VcP) 

1 - ^ 

+ (^ac^M + Sad^bc - Sabred) ( VeTT.) ( VrfVTfe) + (1 - a) J] i^/e^^ + ■ ■ ■ (106) 

Applying incompressible scaling, we are left with the dominant term 



Q, = p + 0{e^). (107) 



1 

+- 




a 




2 




1 




a 


+- 


2 




1 

1 




a 


' 2 



where we have adopted incompressible scaling in the final step. 
The TT component may be written 

Qn. = +5^/({})p + 5^/H{})7r,-5,/(W)V,p-S^/H{c})V.7r, 
+ ls^/{{b, c})V,V.p + Is^J'iic ci})VeV,7r, 

- (1 - a) T,/({6})V,p - (1 - a) T,/HW)Ve7rb 
T^/iib, c}) VbV,p + (2) VVtt, 

U^/'i{b,c}) (V,p) (V,p) + l^t/^/-({c,4) (V,p) (V,vr,) 

f/^/^^({c,d}) (V.vr,) (V.p) + i^f/^/^-({d,e}) (V.tt,) (V.ttJ 

b 

+ (l-a)5^i^,/ef) + --- (108) 

i 

Using the results for the D2Q6 model compiled in Tables [H [2] and [31 this becomes 

Q-Ta = +Sab'n'b - ^^ab^bP + o {^ab^cd + Sac^bd + SadSbc) VcVrfTTfo 

Z o 

- ^ ^ (^^'^ab - 27ra7rb) VfoP - ^ ^ (^afeVTc + (JfocVTa - ^acTTfo) VcTTfo 

6 

+ (l-a)^JVef + (109) 

j 

Applying incompressible scaling, we are left with the dominant terms 

(110) 



= TT - i Vp + i VV - j (tt . Vtt - i V iTT^ + O (e^) . 
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Table 1: The iS'(A^) tensors for the D2Q6 hydrodynamic lattice BGK model 



N = 



5./({})=0 



N 



S/i{a}) = 
Sj\{b}) = 15,, 



N 



N 



S/{{a,b}) = l6ab 
V"({6,c}) = 
S^/{{b,c}) = 

5'^/''({c, d}) = I {dabbed + ^achd + ^adhc) 



N = 3 S/{{a,b,c}) =0 

S/-{{b, C, d}) = i + 6aAd + ^ad^bc) 

S'^/({6, C, d}) = I (5afe4d + SacSbd + SadSbc) 

S^/^i{c,d,e}) = 



S/{{a, b, C, d}) = I {6abScd + ^ac^M + SadSbc) 

S/'^{{b,c,d,e}) = 
S^/i{b,c,d,e}) = 
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Table 2: The T{N) tensors for the D2Q6 hydrodynamic lattice BGK model 

N = T/({}) = 

T/«({}) = 

r./({}) = o 

N=l T/({a}) = 

T/"(W) = 

iV = 2 T/({a, 6}) = ^ (tt^^,, - 27r,7rfe) 

T./({^c})=0 
T,/H{c,4) = 

iV = 3 r/({a,6}) = 

T/''({6,c})=0 
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Table 3: The U{N) tensors for the D2Q6 hydrodynamic lattice BGK model 



N = 





U/P{{}) = 
[/^p-({}) = f/^-p({}) = 
[/^--.({}) = 

K/'''{{}) = u^/^p{{}) = o 

f^TTa [if) - U 




N = 


1 


U/^i{a}) = 
f/^p-({6})=f/^-p({6})=0 
[//''^"({c}) = 

TT '^b'^cfSrlW — l/'AA _lAA AA^ 
'-^■Ka J ~ 2p ['^abOcd + OacObd — OadObc) 


^acT^b) 


N = 


2 


U/-^{{b, c}) = U^^^^P{{b, c}) = {6ab7r, + SacTib - 

t//"'^''({c, d}) = ^ {SacSbd + SadSbc - SabScd) 

K/''i{b,c}) = 
KJ''-^{{d,e}) = 


- SbcT^a) 


N = 


3 


U/P{{a,b,c}) = 
U/^'^{{b,c,d}) = U^^^P{{b,c,d}) =0 
f//''-''({c,rf,e}) =0 
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